{
 "metadata": {},
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rebinning\n",
      "=========\n",
      "\n",
      "Examples showing how to rebin data to produce a smaller or bigger array\n",
      "without (and with) using interpolation.\n",
      "\n",
      "Example 1\n",
      "---------\n",
      "\n",
      "Here we deal with the simplest case where any desired new shape is valid\n",
      "and no interpolation is done on the data to determine the new values. \\*\n",
      "First, floating slices objects are created for each dimension. \\*\n",
      "Second, the coordinates of the new bins are computed from the slices\n",
      "using mgrid. \\* Then, coordinates are transformed to integer indices. \\*\n",
      "And, finally, 'fancy indexing' is used to evaluate the original array at\n",
      "the desired indices."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rebin( a, newshape ):\n",
      "        '''Rebin an array to a new shape.\n",
      "        '''\n",
      "        assert len(a.shape) == len(newshape)\n",
      "\n",
      "        slices = [ slice(0,old, float(old)/new) for old,new in zip(a.shape,newshape) ]\n",
      "        coordinates = mgrid[slices]\n",
      "        indices = coordinates.astype('i')   #choose the biggest smaller integer index\n",
      "        return a[tuple(indices)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we were only interested in reducing the sizes by some integer factor\n",
      "then we could use:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rebin_factor( a, newshape ):\n",
      "        '''Rebin an array to a new shape.\n",
      "        newshape must be a factor of a.shape.\n",
      "        '''\n",
      "        assert len(a.shape) == len(newshape)\n",
      "        assert not sometrue(mod( a.shape, newshape ))\n",
      "\n",
      "        slices = [ slice(None,None, old/new) for old,new in zip(a.shape,newshape) ]\n",
      "        return a[slices]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example 2\n",
      "---------\n",
      "\n",
      "Here is an other way to deal with the reducing case for ndarrays. This\n",
      "acts identically to IDL's rebin command where all values in the original\n",
      "array are summed and divided amongst the entries in the new array. As in\n",
      "IDL, the new shape must be a factor of the old one. The ugly 'evList\n",
      "trick' builds and executes a python command of the form\n",
      "\n",
      "a.reshape(args[0],factor[0],).sum(1)/factor[0]\n",
      "a.reshape(args[0],factor[0],args[1],factor[1],).sum(1).sum(2)/factor[0]/factor[1]\n",
      "\n",
      "etc. This general form is extended to cover the number of required\n",
      "dimensions."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rebin(a, *args):\n",
      "    '''rebin ndarray data into a smaller ndarray of the same rank whose dimensions\n",
      "    are factors of the original dimensions. eg. An array with 6 columns and 4 rows\n",
      "    can be reduced to have 6,3,2 or 1 columns and 4,2 or 1 rows.\n",
      "    example usages:\n",
      "    >>> a=rand(6,4); b=rebin(a,3,2)\n",
      "    >>> a=rand(6); b=rebin(a,2)\n",
      "    '''\n",
      "    shape = a.shape\n",
      "    lenShape = len(shape)\n",
      "    factor = asarray(shape)/asarray(args)\n",
      "    evList = ['a.reshape('] + \\\n",
      "             ['args[%d],factor[%d],'%(i,i) for i in range(lenShape)] + \\\n",
      "             [')'] + ['.sum(%d)'%(i+1) for i in range(lenShape)] + \\\n",
      "             ['/factor[%d]'%i for i in range(lenShape)]\n",
      "    print ''.join(evList)\n",
      "    return eval(''.join(evList))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The above code returns an array of the same type as the input array. If\n",
      "the input is an integer array, the output values will be rounded down.\n",
      "If you want a float array which correctly averages the input values\n",
      "without rounding, you can do the following instead.\n",
      "\n",
      "a.reshape(args[0],factor[0],).mean(1) <BR>\n",
      "a.reshape(args[0],factor[0],args[1],factor[1],).mean(1).mean(2)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rebin(a, *args):\n",
      "    shape = a.shape\n",
      "    lenShape = len(shape)\n",
      "    factor = asarray(shape)/asarray(args)\n",
      "    evList = ['a.reshape('] + \\\n",
      "             ['args[%d],factor[%d],'%(i,i) for i in range(lenShape)] + \\\n",
      "             [')'] + ['.mean(%d)'%(i+1) for i in range(lenShape)]\n",
      "    print ''.join(evList)\n",
      "    return eval(''.join(evList))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Some test cases:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 2-D case\n",
      "a=rand(6,4)\n",
      "print a\n",
      "b=rebin(a,6,4)\n",
      "print b\n",
      "b=rebin(a,6,2)\n",
      "print b\n",
      "b=rebin(a,3,2)\n",
      "print b\n",
      "b=rebin(a,1,1)\n",
      "\n",
      "# 1-D case\n",
      "print b\n",
      "a=rand(4)\n",
      "print a\n",
      "b=rebin(a,4)\n",
      "print b\n",
      "b=rebin(a,2)\n",
      "print b\n",
      "b=rebin(a,1)\n",
      "print b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example 3\n",
      "---------\n",
      "\n",
      "A python version of congrid, used in IDL, for resampling of data to\n",
      "arbitrary sizes, using a variety of nearest-neighbour and interpolation\n",
      "routines."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as n\n",
      "import scipy.interpolate\n",
      "import scipy.ndimage\n",
      "\n",
      "def congrid(a, newdims, method='linear', centre=False, minusone=False):\n",
      "    '''Arbitrary resampling of source array to new dimension sizes.\n",
      "    Currently only supports maintaining the same number of dimensions.\n",
      "    To use 1-D arrays, first promote them to shape (x,1).\n",
      "    \n",
      "    Uses the same parameters and creates the same co-ordinate lookup points\n",
      "    as IDL''s congrid routine, which apparently originally came from a VAX/VMS\n",
      "    routine of the same name.\n",
      "\n",
      "    method:\n",
      "    neighbour - closest value from original data\n",
      "    nearest and linear - uses n x 1-D interpolations using\n",
      "                         scipy.interpolate.interp1d\n",
      "    (see Numerical Recipes for validity of use of n 1-D interpolations)\n",
      "    spline - uses ndimage.map_coordinates\n",
      "\n",
      "    centre:\n",
      "    True - interpolation points are at the centres of the bins\n",
      "    False - points are at the front edge of the bin\n",
      "\n",
      "    minusone:\n",
      "    For example- inarray.shape = (i,j) & new dimensions = (x,y)\n",
      "    False - inarray is resampled by factors of (i/x) * (j/y)\n",
      "    True - inarray is resampled by(i-1)/(x-1) * (j-1)/(y-1)\n",
      "    This prevents extrapolation one element beyond bounds of input array.\n",
      "    '''\n",
      "    if not a.dtype in [n.float64, n.float32]:\n",
      "        a = n.cast[float](a)\n",
      "    \n",
      "    m1 = n.cast[int](minusone)\n",
      "    ofs = n.cast[int](centre) * 0.5\n",
      "    old = n.array( a.shape )\n",
      "    ndims = len( a.shape )\n",
      "    if len( newdims ) != ndims:\n",
      "        print \"[congrid] dimensions error. \" \\\n",
      "              \"This routine currently only support \" \\\n",
      "              \"rebinning to the same number of dimensions.\"\n",
      "        return None\n",
      "    newdims = n.asarray( newdims, dtype=float )    \n",
      "    dimlist = []\n",
      "\n",
      "    if method == 'neighbour':\n",
      "        for i in range( ndims ):\n",
      "            base = n.indices(newdims)[i]\n",
      "            dimlist.append( (old[i] - m1) / (newdims[i] - m1) \\\n",
      "                            * (base + ofs) - ofs )\n",
      "        cd = n.array( dimlist ).round().astype(int)\n",
      "        newa = a[list( cd )]\n",
      "        return newa\n",
      "    \n",
      "    elif method in ['nearest','linear']:\n",
      "        # calculate new dims\n",
      "        for i in range( ndims ):\n",
      "            base = n.arange( newdims[i] )\n",
      "            dimlist.append( (old[i] - m1) / (newdims[i] - m1) \\\n",
      "                            * (base + ofs) - ofs )\n",
      "        # specify old dims\n",
      "        olddims = [n.arange(i, dtype = n.float) for i in list( a.shape )]\n",
      "\n",
      "        # first interpolation - for ndims = any\n",
      "        mint = scipy.interpolate.interp1d( olddims[-1], a, kind=method )\n",
      "        newa = mint( dimlist[-1] )\n",
      "\n",
      "        trorder = [ndims - 1] + range( ndims - 1 )\n",
      "        for i in range( ndims - 2, -1, -1 ):\n",
      "            newa = newa.transpose( trorder )\n",
      "\n",
      "            mint = scipy.interpolate.interp1d( olddims[i], newa, kind=method )\n",
      "            newa = mint( dimlist[i] )\n",
      "\n",
      "        if ndims > 1:\n",
      "            # need one more transpose to return to original dimensions\n",
      "            newa = newa.transpose( trorder )\n",
      "\n",
      "        return newa\n",
      "    elif method in ['spline']:\n",
      "        oslices = [ slice(0,j) for j in old ]\n",
      "        oldcoords = n.ogrid[oslices]\n",
      "        nslices = [ slice(0,j) for j in list(newdims) ]\n",
      "        newcoords = n.mgrid[nslices]\n",
      "\n",
      "        newcoords_dims = range(n.rank(newcoords))\n",
      "        #make first index last\n",
      "        newcoords_dims.append(newcoords_dims.pop(0))\n",
      "        newcoords_tr = newcoords.transpose(newcoords_dims)\n",
      "        # makes a view that affects newcoords\n",
      "\n",
      "        newcoords_tr += ofs        \n",
      "\n",
      "        deltas = (n.asarray(old) - m1) / (newdims - m1)\n",
      "        newcoords_tr *= deltas\n",
      "\n",
      "        newcoords_tr -= ofs\n",
      "\n",
      "        newa = scipy.ndimage.map_coordinates(a, newcoords)\n",
      "        return newa\n",
      "    else:\n",
      "        print \"Congrid error: Unrecognized interpolation type.\\n\", \\\n",
      "              \"Currently only \\'neighbour\\', \\'nearest\\',\\'linear\\',\", \\\n",
      "              \"and \\'spline\\' are supported.\"\n",
      "        return None"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}